\section{Ice Thermal Storage }\label{ice-thermal-storage}

\subsection{Simple Ice Storage Model}\label{simple-ice-storage-model}

\subsubsection{Overview}\label{overview-020}

This model (object name ThermalStorage:Ice:Simple) simulates the performance of an ice storage tank. The model is based on an integrated chiller and storage tank model developed for a special optimization project (Henze and Krarti 2002). This implementation removed the integrated chiller and allows the tank to be charged by any chiller configuration in an EnergyPlus plant loop..

This thermal storage model is based on a simple simulation of an ice storage tank with a fixed capacity.~ The tank is charged, or frozen, in an ice-on-coil configuration where ice builds up on the outside of the tubes carrying the brine or glycol solution from the chiller.~ There are two discharge (melt) options, internal or external.~ Internal melt uses the same fluid tubes for charging and discharging.~ External melt uses a separate fluid path for discharge such that the outer layers of ice melt first.~ The ice storage model includes an implied 3-way valve to control the amount if charge/discharge based on the incoming water temperature and the outlet node setpoint temperature.~ The tank is controlled using the PlantEquipmentOperation:ComponentSetpoint plant operation scheme, and requires that a setpoint be placed by a setpoint manager on the ice storage Plant Outlet Node.~ The model currently does not calculate any tank losses.

\subsubsection{Model Description}\label{model-description-013}

\paragraph{Initialization}\label{initialization-000}

The storage tank is assumed to be fully charged (full of ice) at the beginning of each environment.~ The tank is then allowed to charge and discharge during the warmup days of the environment.

\paragraph{Control}\label{control-000}

The ice storage module is a passive component and will accept any flow rate placed on the inlet node.~ The load on the ice storage tank is determined by the inlet node flow rate, the inlet node temperature, and the outlet node setpoint temperature.~ A positive load indicates a request for cooling and the tank discharges if possible.~ A negative load indicates a request for charging the tank if possible.~ A zero load indicates the tank is dormant in which case all flow bypasses the tank and the outlet node temperature is set to the inlet node temperature.~ The model includes an implied 3-way valve which controls the leaving water temperature to the outlet node setpoint temperature if possible.

\paragraph{Charging}\label{charging}

When charging is requested, the following limits are calculated to determine the actual charging rate:

\begin{enumerate}
\item If the entering water temperature is greater than -1\(^{\circ}\)C, the charging rate is zero.
\item If the entering water temperature is greater than or equal to the outlet setpoint temperature, the charging rate is zero.
\item If the current fraction of ice stored is 1.0, the tank is fully charge, and the charging rate is zero.
\item If the requested charging rate times the current system timestep is more than the remaining uncharged storage capacity, the charging rate is limited to that required to completely fill the tank.
\item The maximum charging rate which the chilled water flow can provide is determined by the entering water temperature and flow rate and an outlet temperature which is the minimum of -1\(^{\circ}\)C or the outlet setpoint temperature.
\item The maximum charging rate which the tank can accept is calculated by the following equations developed in the prior work (Henze and Krarti 2002):
\end{enumerate}

\begin{equation}
\begin{split}
{UAIceCh} = (& 1.3879 - 7.6333*{y} + 26.3423*{y^2} - 47.6084*{y^3} + 41.8498*{y^4} \\
&- 14.2948*{y^5} ) * {ITSNomCap} / {TimeInterval} / 10.0
\end{split}
\end{equation}

where:

\emph{UAIceCh} is the UA value for charging

\emph{y} is the current ice fraction stored

\emph{ITSNomCap} is the nominal storage capacity (GJ)

TimeInterval is 3600 (s).

The smallest charging rate determined by the above rules is selected and the corresponding leaving water temperature is calculated.

\paragraph{Discharging}\label{discharging}

When discharging is requested, the following limits are calculated to determine the actual charging rate:

\begin{enumerate}
\item If the entering water temperature is less than 1\(^{\circ}\)C, the discharge rate is zero.
\item If the entering water temperature is less than or equal to the outlet setpoint temperature, the discharge rate is zero.
\item If the current fraction of ice stored is 0.0, the tank is fully discharged, and the discharge rate is zero.
\item If the requested discharge rate times the current system timestep is more than the remaining charged storage capacity, the discharge rate is limited to that required to completely deplete the ice in the tank.
\item The maximum discharge rate which the chilled water flow can accept is determined by the entering water temperature and flow rate and an outlet temperature which is the maximum of 1\(^{\circ}\)C or the outlet setpoint temperature.
\item The maximum discharge rate which the tank can deliver is calculated by the following equations developed in the prior work (Henze and Krarti 2002):
\end{enumerate}

For ice-on-coil internal melt:

\begin{equation}
\begin{split}
{UAIceDisCh} = (& 1.3879 - 7.6333*{y} + 26.3423*{y^2} - 47.6084*{y^3} + 41.8498*{y^4} \\
&- 14.2948*{y^5} ) * {ITSNomCap} / {TimeInterval} / 10.0
\end{split}
\end{equation}

For ice-on-coil external melt:

\begin{equation}
\begin{split}
{UAIceDisCh} = (& 1.1756 - 5.3689*{y} + 17.3602*{y^2} - 30.1077*{y^3} + 25.6387*{y^4} \\
&- 8.5102*{y^5} ) * {ITSNomCap} / {TimeInterval} / 10.0
\end{split}
\end{equation}

where:

UAIceDisCh is the UA value for discharging

\emph{y} = 1 - current ice fraction stored

\emph{ITSNomCap} is the nominal storage capacity (GJ)

TimeInterval is 3600 (s).

The smallest discharge rate determined by the above rules is selected and the corresponding leaving water temperature is calculated.

\subsubsection{References}\label{references-028}

Henze, Gregor P. and Moncef Krarti. 2002. Predictive Optimal Control of Active and Passive Building Thermal Storage Inventory, Final Report for Phase I:~ Analysis, Modeling, and Simulation.~ U.S. Department of Energy National Energy Technology Laboratory Cooperative Agreement DE-FC-26-01NT41255, December 2002.

\subsection{Detailed Ice Storage Model}\label{detailed-ice-storage-model}

The following section describes how the detailed ice storage model works in EnergyPlus (object name ThermalStorage:Ice:Detailed).

\subsubsection{Charging and Discharging Equation}\label{charging-and-discharging-equation}

The actual performance of the ice storage unit depends on the physical geometry, materials, and characteristics of the ice storage unit. Modeling both the charging and discharging performance of a particular ice storage unit is accomplished in EnergyPlus using Curve objects that establish the relationships between various parameters of the unit and the output of the unit itself. The input that controls what curve objects are used and what parameters are the independent variables in those curves are: Discharging Curve Variable Specifications, Discharging Curve Name, Charging Curve Variable Specifications, and Charging Curve Name.

The Variable Specifications input fields determine which parameters are used as the independent variables in the curves referenced by the Curve Name inputs.  There are four options for the Variable Specifications fields: FractionChargedLMTD, FractionDischargedLMTD, LMTDMassFlow, and LMTDFractionCharged.  Each option defines which two variables are the independent variables in the curve equation as well as the order.  For example, FractionChargedLMTD uses the fraction of full charge of the storage component as the first variable and the normalized LMTD (or LMTD*) as the second variable.  Note that using options that use either LMTD or MassFlow actually use normalized values LMTD* and MassFlow*.  LMTD is normalized by dividing the actual value by 10\(^{\circ}\)C while MassFlow is normalized by dividing the actual mass flow rate by the equivalent of 100 gallons per minute.

One equation that might be used to characterize the discharging performance of an ice storage component is a QuadraticLinear Curve.  If the discharging equation referenced in the detailed ice storage unit is a QuadraticLinear equation that uses FractionDischargedLMTD for the variable specification, then this would result in the following equation for characterizing the ice storage component's discharging process:

\begin{equation}
{q^*} = \left[ {C1 + C2\left( {1 - {P_c}} \right) + C3{{\left( {1 - {P_c}} \right)}^2}} \right] + \left[ {C4 + C5\left( {1 - {P_c}} \right) + C6{{\left( {1 - {P_c}} \right)}^2}} \right]\Delta T_{lm}^*
\end{equation}

or

\begin{equation}
{q^*} = \left[ {C1 + C2\left( {{P_d}} \right) + C3{{\left( {{P_d}} \right)}^2}} \right] + \left[ {C4 + C5\left( {{P_d}} \right) + C6{{\left( {{P_d}} \right)}^2}} \right]\Delta T_{lm}^*
\end{equation}

where:

\begin{equation}
{q^*} \equiv \frac{{q\Delta t}}{{{Q_{stor}}}}
\end{equation}

\begin{equation}
\Delta T_{lm}^* \equiv \frac{{\Delta {T_{lm}}}}{{\Delta {T_{nominal}}}}
\end{equation}

\begin{equation}
\Delta {T_{lm}} \equiv \frac{{{T_{brine,in}} - {T_{brine,out}}}}{{\ln \left( {\frac{{{T_{brine,in}} - {T_{brine,freeze}}}}{{{T_{brine,out}} - {T_{brine,freeze}}}}} \right)}}
\end{equation}

q is the instantaneous heat transfer rate

Q\(_{stor}\) is the total latent storage capacity

\(\Delta\)t is a time step used in the curve fit (usually one hour)

\(\Delta\)T\(_{nominal}\) is a nominal temperature difference (18\(^{\circ}\)F = 10\(^{\circ}\)C)

T\(_{brine,in}\) is the tank brine inlet temperature

T\(_{brine,out}\) is the tank brine outlet temperature

T\(_{freeze}\) is the freezing temperature of water or the latent energy storage material

P\(_{c}\) is the fraction charged

P\(_{d}\) is the fraction discharged.

Likewise, if the charging equation referenced in the detailed ice storage unit is a QuadraticLinear equation that uses FractionChargedLMTD for the variable specification, then this would result in the following equation for characterizing the ice storage component's charging process:

\begin{equation}
{q^*} = \left[ {C1 + C2\left( {{P_c}} \right) + C3{{\left( {{P_c}} \right)}^2}} \right] + \left[ {C4 + C5\left( {{P_c}} \right) + C6{{\left( {{P_c}} \right)}^2}} \right]\Delta T_{lm}^*
\end{equation}

Note that the time step used for the curve fit of performance data might differ from the time step used within the EnergyPlus simulation.~ These are actually two separate time steps and are kept separate.

\subsubsection{Charging Algorithm}\label{charging-algorithm}

During charging, manufacturers have stated that they attempt to charge the unit at the maximum rate until the unit is completely charged.~ This, of course, occurs during off-peak electric hours.~ Thus, once the setpoint has been scheduled for charging, the unit will charge at the maximum possible rate.~ This means that the flow rate through the ice storage device equals the flow to the component (or no bypass).~ The only time flow to the ice storage unit would be reduced is at the end of the charge cycle when more ice making capacity is available in a particular time step than is needed to fully charge the tank.~ In this case, the flow to the tank would be reduced appropriately to top off the tank storage capacity. ~We also have a setpoint goal for the outlet temperature of the ice storage device as defined by the setpoint schedule.

In solving the performance of the ice storage unit, we have effectively two equations.  One of these equations is determined by the user input and characterizes the output of the unit (\(q^*\)) as a function of LMTD* and fraction charged/discharged or normalized mass flow rate. The other equation is:

\begin{equation}
q = \dot m{C_p}\left( {{T_i} - {T_o}} \right)
\end{equation}

Both of these equations have q and T\(_{o}\) as unknowns.~ However, since the setpoint temperature is the goal for T\(_{o}\), we can use this as an initial guess for T\(_{o}\).~ Below is an outline of the algorithm:

\begin{itemize}
  \item
    Initialize T\(_{o}\) = T\(_{set}\)
  \item
    Calculate LMTD*
  \item
    Calculate q* from charging equation for the current percent charged (We will assume that the EnergyPlus time step is sufficiently small so that we do not need to find the average percent charged for the time step. This was necessary when one hour time steps were used as in BLAST, but EnergyPlus generally uses relatively short time steps. Since there is already some iteration involved in the solution, we would like to avoid another layer of iteration if at all possible. One alternative that could be implemented would be to make a second pass with a closer average value based on what happens during the time step. This would effectively double the execution time for the model and would need to be justified before implementation.)
  \item
    Calculate T\(_{o,new}\) and compare it to T\(_{o}\)
  \item
    Use T\(_{o,new}\) to calculate a new LMTD* and iterate until T\(_{o}\) converges to some acceptable level
\end{itemize}

Charging would continue in subsequent time steps until the final state of the ice storage unit at the end of a particular time step is fully charged. If running a chiller would ``overcharge'' the tank, then the flow to the tank would be reduced (greater than zero bypass flow) while maintaining the same setpoint temperature coming out of the tank (though not necessarily out of the component).

\subsubsection{Discharging Algorithm}\label{discharging-algorithm}

During discharging, we cannot assume that all of the flow is sent through the ice storage unit and thus some of it may be bypassed around it locally. This ice storage model includes a built-in bypass leg to accommodate this without requiring the user to enter this additional information. This also allows the bypass leg/valve to be controlled by the ice storage unit.

While we cannot assume that all of the flow is sent through the ice storage unit, we can use that as an initial guess in order to determine the current performance of the ice storage system. Most of the discharging algorithm then becomes very similar to the charging process.

In solving the performance of the ice storage unit, we have effectively two equations.  One of these equations is determined by the user input and characterizes the output of the unit (\(q^*\)) as a function of LMTD* and fraction charged/discharged or normalized mass flow rate. The other equation is:

\begin{equation}
q = \dot m{C_p}\left( {{T_i} - {T_o}} \right)
\end{equation}

Both of these equations have q and T\(_{o}\) as unknowns.~ However, since the setpoint temperature is the goal for T\(_{o}\), we can use this as an initial guess for T\(_{o}\).~ Below is an outline of the algorithm:

\begin{itemize}
  \item
    Initialize T\(_{o}\) = T\(_{set}\)
  \item
    Calculate LMTD*
  \item
    Calculate q* from charging equation for the current percent charged (we will assume that the EnergyPlus time step is sufficiently small so that we do not need to find the average percent charged for the time step; this was necessary when one hour time steps were used as in BLAST, but EnergyPlus generally uses relatively short time steps)
  \item
    Calculate T\(_{o,new}\) and compare it to T\(_{o}\)
  \item
    Use T\(_{o,new}\) to calculate a new LMTD* and iterate until T\(_{o}\) converges to some acceptable level
  \item
    Once T\(_{o}\) has converged, we need to compare this value again to T\(_{set}\).~ If T\(_{o}\) is greater than or equal to T\(_{set}\), then we either just met the load (equal) or can't quite meet it (greater than).~ In this case, we simply leave T\(_{o}\) as is because we cannot meet the setpoint temperature.
  \item
    If T\(_{o}\) is less than T\(_{set}\), then we have more capacity available than we need.~ In this case, we need to bypass some of the flow.~ Since the load on the ice storage device and the outlet temperature are not changing (we are just reducing the flow), we only need to split the flow and do not need to recalculate the action of the ice storage device.~ Some systems may be slightly dependent on the actual flow through the device itself.~ However, in an actual application, this only means that a slightly different amount will bypass the device.~ The overall energy impact will be the same and thus it is not necessary to be concerned about flow rate dependence.
\end{itemize}

Discharging would continue in subsequent time steps until the final state of the ice storage unit at the end of a particular time step is fully discharged.

\subsubsection{References}\label{references-1-011}

Strand, R.K. 1992. ``Indirect Ice Storage System Simulation,'' M.S. Thesis, Department of Mechanical and Industrial Engineering, University of Illinois at Urbana-Champaign.

\subsection{Other Ice Thermal Storage Options}

Thermal storage can also be modeled as integrated with a packaged air-conditioning unit. For this option, see Section \ref{packaged-thermal-storage-cooling-coil}.

\subsection{PCM Thermal Storage Model}\label{pcm-thermal-storage-model}

\subsubsection{Overview}\label{pcm-overview}

This model (object name \texttt{ThermalStorage:PCM}) represents a well-mixed phase-change material (PCM) tank connected to a plant loop on two sides: a \emph{plant side} (used to charge the PCM) and a \emph{use side} (used to discharge the PCM to meet chilled-water load). The tank stores and releases \emph{latent} energy as the PCM melts and freezes over a finite temperature range around the material’s peak melting temperature. The model is intentionally simple and robust for whole-building simulations and uses effectiveness-based heat-exchanger relationships for outlet temperatures together with an energy balance to track the tank’s state of charge (SOC).

Key features and assumptions:

\begin{itemize}
  \item The tank is \emph{fully mixed} (no stratification), so one bulk PCM temperature represents the entire volume.
  \item Only \emph{one side} (plant or use) is active in a timestep, selected by a hysteresis SOC controller; simultaneous charge and discharge are not modeled.
  \item Heat loss to ambient is a constant rate input (W) applied every system timestep.
  \item PCM thermophysical behavior is captured by a material function that maps temperature to specific enthalpy and vice versa; the phase change occurs over a finite temperature interval around the peak melting temperature.
  \item The tank’s energy state is tracked in joules (J). Heat-transfer rates are in watts (W).
\end{itemize}

\subsubsection{Model Description}\label{pcm-model-description}

\paragraph{Variables and parameters}\label{pcm-variables}

\begin{tabular}{ll}
$E$ & Stored energy in the PCM tank (J) \\
$E_{\max} = m_{\mathrm{PCM}} \, h_{\ell}$ & Maximum latent storage (J), mass $\times$ latent heat \\
$m_{\mathrm{PCM}}$ & Tank capacity as PCM mass (kg) \\
$h_{\ell}$ & Latent heat of fusion (J/kg) \\
$T_\mathrm{m}$ & Peak melting temperature of the PCM (°C) \\
$\Delta T_{\ell}, \Delta T_{h}$ & Lower/upper half-widths of the melting range (°C) \\
$\epsilon$ & Heat-exchanger effectiveness (–) \\
$\dot m_\mathrm{use}, \dot m_\mathrm{pl}$ & Mass flow rates on use/plant sides (kg/s) \\
$c_p$ & Water specific heat (J/(kg·K)) \\
$\dot Q_\mathrm{loss}$ & Constant heat loss to ambient (W) \\
$\mathrm{SOC} = E/E_{\max}$ & State of charge (–) \\
$\mathrm{SOC}_{\mathrm{in}}, \mathrm{SOC}_{\mathrm{out}}$ & Hysteresis thresholds for charge start/stop (–)
\end{tabular}

\paragraph{Initialization}\label{pcm-initialization}

At the beginning of each environment, the tank may be initialized to a specified initial charge fraction (default $=1.0$), i.e.,
\begin{equation}
E(t_0) = \mathrm{SOC}_0 \, E_{\max}.
\end{equation}
If warmup days are used, the tank can charge/discharge during warmup; thus the SOC at the first reported timestep may be lower than the initial value unless the model is explicitly reset when warmup ends.

\paragraph{Outlet temperature calculations}\label{pcm-outlet-temps}

For robustness, outlet water temperatures are computed with an effectiveness relation relative to the PCM phase-change temperature band. Let $T_{\mathrm{in,pl}}$ and $T_{\mathrm{in,use}}$ be inlet temperatures. The model forms target outlet temperatures
\begin{align}
T_{\mathrm{out,pl}} &= T_{\mathrm{in,pl}} - \epsilon \,\bigl(T_{\mathrm{in,pl}} - T_f\bigr), \\
T_{\mathrm{out,use}} &= T_{\mathrm{in,use}} + \epsilon \,\bigl(T_m - T_{\mathrm{in,use}}\bigr),
\end{align}
where $T_f$ and $T_m$ are representative freezing and melting temperatures (derived from the PCM material data around $T_\mathrm{m}$). These targets are then clipped to plant loop temperature bounds. The corresponding side heat-transfer rates are
\begin{align}
\dot Q_{\mathrm{pl}} &= \dot m_{\mathrm{pl}} \, c_p \, \bigl(T_{\mathrm{in,pl}} - T_{\mathrm{out,pl}}\bigr), \\
\dot Q_{\mathrm{use}} &= \dot m_{\mathrm{use}} \, c_p \, \bigl(T_{\mathrm{in,use}} - T_{\mathrm{out,use}}\bigr).
\end{align}

\paragraph{SOC control and flow selection}\label{pcm-control}

A simple hysteresis controller determines which side is active:
\begin{itemize}
  \item If $\mathrm{SOC} \le \mathrm{SOC}_{\mathrm{in}}$ (e.g., $0.40$), \emph{charging mode} is enabled and the component requests only plant-side flow: $\dot m_{\mathrm{pl}} = \dot m_{\mathrm{pl,design}}$, $\dot m_{\mathrm{use}} = 0$.
  \item If $\mathrm{SOC} \ge \mathrm{SOC}_{\mathrm{out}}$ (e.g., $0.95$), \emph{discharging mode} is enabled and the component requests only use-side flow: $\dot m_{\mathrm{use}} = \dot m_{\mathrm{use,design}}$, $\dot m_{\mathrm{pl}} = 0$.
\end{itemize}
Only one side is active per timestep by construction. If the plant loop solver does not honor a requested flow (e.g., pump off), the actual heat exchange on that side is zero for that step.

\paragraph{Energy update}\label{pcm-energy-update}

Let $\Delta t$ be the system timestep in seconds. The net power into the PCM is
\begin{equation}
\dot Q_{\mathrm{net}} = \dot Q_{\mathrm{pl}} + \dot Q_{\mathrm{use}} - \dot Q_{\mathrm{loss}}.
\end{equation}
The stored energy state is advanced explicitly:
\begin{equation}
E(t+\Delta t) = \min\!\Bigl(E_{\max}, \max\!\bigl(0,\, E(t) + \dot Q_{\mathrm{net}} \,\Delta t \bigr)\Bigr),
\end{equation}
and $\mathrm{SOC}(t+\Delta t) = E(t+\Delta t)/E_{\max}$.

\paragraph{PCM temperature from enthalpy}\label{pcm-temp-from-enthalpy}

The bulk PCM temperature is diagnosed from $E$ via the specific enthalpy function of the PCM material, $h(T)$, which includes both sensible and latent contributions across the melting range:
\begin{equation}
h(T) = h_\mathrm{ref} + \int_{T_\mathrm{ref}}^{T} c_{p,\mathrm{PCM}}(\theta)\, d\theta + \Phi(\theta; T_\mathrm{m}, \Delta T_{\ell}, \Delta T_h) \, h_{\ell}.
\end{equation}
Here $\Phi\in[0,1]$ is a smooth phase-fraction function centered at $T_\mathrm{m}$ with the specified half-widths. The model inverts $h(T)$ numerically (bisection) to obtain $T$ from the target mass-specific enthalpy $E/m_{\mathrm{PCM}}$.

\paragraph{Limits and notes}\label{pcm-limits}

\begin{itemize}
  \item The tank is treated as fully mixed; axial stratification is not modeled.
  \item Only one side (charge or discharge) may be active per timestep.
  \item The effectiveness model is empirical; it does not resolve detailed HX geometry.
  \item A constant heat-loss rate is assumed; radiative/convective breakdown is not modeled.
  \item Accurate PCM material data are required for meaningful phase-transition behavior.
\end{itemize}

\subsubsection{Performance Checks}\label{pcm-performance-checks}

For each system timestep, the following diagnostic identity should hold within numerical tolerance:
\begin{equation}
\Delta E \;\approx\; \bigl(\dot Q_{\mathrm{pl}} + \dot Q_{\mathrm{use}} - \dot Q_{\mathrm{loss}}\bigr)\,\Delta t.
\end{equation}
Large discrepancies indicate unmet flow requests or inconsistent reporting configurations.

\subsubsection{Inputs and Outputs}\label{pcm-io-refs}

See the I/O Reference section \emph{ThermalStorage:PCM} for field definitions (availability schedule, node names, PCM material, tank capacity, heat loss, and design flow rates). Typical output variables include:
\begin{itemize}
  \item \texttt{Thermal Energy Storage Energy Stored} \, (J)
  \item \texttt{Thermal Energy Storage Percent Capacity} \, (\%)
  \item \texttt{Thermal Energy Storage Plant Side Heat Transfer Rate} \, (W)
  \item \texttt{Thermal Energy Storage Use Side Heat Transfer Rate} \, (W)
  \item \texttt{Thermal Energy Storage Heat Loss Rate} \, (W)
  \item \texttt{Thermal Energy Storage Tank Temperature} \, (°C)
\end{itemize}
